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VESSEL CENTERLINE DETERMINATION 
FIELD OF THE INVENTION 
The present invention relates to the processing of data sets, for example medical image 
data sets obtained from CT images. 
5 BACKGROUND OF THE INVENTION 

CT imaging systems can be used for checking blood vessels, especially for detecting 
narrowing (stenosis) in the blood vessels. Fig, I is a schematic illustration of a display 
showing two blood vessel views, a view 100 showing a blood vessel 102 and a view 122, a 
slice generally perpendicular to vessel 102, along a plane marked as 120. It should be noted 
10 that in most situations, vessel 102 will not lie unifonnly in one plane. It is desirable to 
determine a centerline 108 of vessel 102, so that slices 122 can be made perpendicular thereto 
and allowing various views along and of the blood vessel. For example, the vessel may be 
checked for narrowing by providing a series of slices perpendicular to the blood vessel or the 

■ 

vessel may be spread out into a flat image, for inspecting its wall 
15 An exemplary current practice is to have a user mark multiple points (e,g,, 16-32 

points) along an estimated centerlme 108 and then connect them with straight lines or use 

various algorithmic methods to obtain a more exact determination of the centerline. 

It should be noted that centerline determination is generally not a trivial task, In many 

cases parts of the blood vessel are missing, have non-uniform CT number (Hounsfleld 
20 numbers) and/or the edges of the blood vessels are not clear. Exemplary reasons for this are: 

(a) Even if the blood vessel is imaged using a contrast material, this material may not 
be unifonnly distributed along the blood vessel. 

(b) Partial volume effects, especially near bones. 

(c) Nearby tissue with similar absorption (e.g., cortical bone, mairow and kidney 
25 tissue). For example, a boundary 118 might not be visible next to a bone 116 

(d) Narrowing, splitting and/or other geometrical properties of vessels. 

(e) Some vessels contain stents. 

(f) Nearby vessels may appear to meet and merge and (hen diverge. 

(g) Various effects may cause a vessel to appear to include loops. 

30 (h) An incorrect centerline can cause the appearance of a narrowing while viewing the 

vessel. 

0 An unsmooth centerline can cause difficulty in viewing and/or understanding the 
entire blood vessel 

1 
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(j) Occlusions. 

SUMMARY OF THE INVENTION 

Ad aspect of some embodiments of the invention relates to a method of centerline 
determination, in which a blood vessel volume is identified based on an inaccurate internal 
path and then a more accurate centeriine is determined 

An aspect of some embodiments of the invention relates to a targeted marching 
algorithm based method! in which a cost function of a point depends both on a local cost and 
on an estimated cost to reach a target. In an exemplary application, a contour is expanded from 
one or more starting points to a target point, with a heap or other data structure being used to 
store points and their associated cost Once the contour passes a point, it is marked "alive" and 
taken out of the heap. In an exemplary embodiment of the invention, the estimated cost to 
target for a current point is based on a function of the geometrical distance to target and the 
average cost (per unit length) to (he current point Optionally, this cost is normalized, for 
example to equate between different starting points. 

Optionally, if an "alive" point is reached by the contour a second time with a lower cost 
than before, it may be reinserted into the heap for further consideration. In an exemplary 
embodiment of the invention, the cost function is, at least most of the time, an estimate from 
below the real cost 

In an exemplary embodiment of the invention, targeted marching is used for finding an 
initial possibly inaccurate centeriine (which is actually an internal path which may tend to the 
center). Optionally, one or more of the following properties are desired in the method used for 
finding this internal path, and may be provided by some implementations of targeted marching: 

(a) geometrically correct distance calculations, resulting from a greater consistency 
with underlying continuous data than provided by other methods (e.g., as data is sampled finer, 
the accuracy of measurement increases); in one example correct diagonal distance calculations 
are provided; 

(b) providing a path within the volume of the blood vessel; 

(c) ability to woric with a relatively small nnmber of starting points; and 

+ 

(d) focusing a search propagation to within the blood vessel volume. 

In other embodiments of the invention, a different path finding algorithm is used. 
In some embodiments of the invention, targeted marching is used to find a mom exact 
center line, once the blood vessel volume is identified. 

2 
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An aspect of some embodiments of the invention relates to a method of segmentation 
of a blood vessel to separate out the data points of a blood vessel from a complete data set, in 
which different sections along the blood vessels are segmented from the data set using 
different local histograms for bloodvessel data values. Optionally, the histograms are averaged 

5 for nearby segments. 

An aspect of some embodiments of the invention relates to a segmentation method in 
which a curvature limitation is applied during point propagation. Optionally, this curvature 
limitation trades-off the possibility of missing part of the vessel with the ability to reduce non- 
uniform leakage and/or bridge gaps. 

1 0 An aspect of some embodiments of the invention relates to cost functions which reduce 

escape of a search point (e.g., at the edge of a contour) from a volume of a blood vessel. In an 
exemplary embodiment of the invention, a path finding method in a blood vessel uses a cost 
function for a point based on a distribution of data values, to determine a probability of the 
point being inside or outside of the blood vessel and/or an estimation of curvature. In an 

15 exemplary embodiment of the invention, a cost function includes a factor of an angle 
limitation, reducing angular motion of the search point relative to a starting point. In an 
exemplary embodiment of the invention, a cost function includes a distance limitation based 
on an average distance in the cost metric found in nearby sections of the blood vessaL In an 
exemplary embodiment of the invention, a cost function includes limiting the propagation of 

20 points from a same source that propagate too for (e.g. a based on a parameter value). In an 
exemplary embodiment of the invention, a cost Amotion includes a limitation that prevent 
forward (along the blood vessel) propagation, which also tends to limit leakage. Optionally, 
the leakage is reduced such that fewer than 50% fewer than 20% or fewer than 10% of the 
points considered by the segmentation method are leakage points, for a typical execution, for 

25 example on 50% of images tor which a centerline is sought. 

An aspect of some embodiments of the invention relates to propagation of 
parameterization in which the parameter values are propagated in a manner which is uniform, 
smooth and/or generally parallel to a starting line. 

An aspect of some embodiments of the invention relates to a method of finding a 

30 centerline, in a manner which trades off curvature (of the centerline) and staying in a center. 
Optionally, the method is generally indifferent to variations in vessel diameter and/or gives 
good results for large diameter vessels. 
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An aspect of some embodiments of the invention relates to centerline and/or 
segmentation finding method which has a typically CPU usage of 0(nlogn) and/or a typical 
memory usage of O(n), where n is the number of points actually found inside the blood vessel 
section of interest. In an exemplary embodiment of the invention, these limits are available due 
5 to restricting leakage outside the blood vessel. In some cases an additional fixed cost is added, 
for example o(N)» where N is the number of points in the whole data set, for forming of a 
histogram. Alternatively, a statistical histogram may be created, for example using only the 
same number of points as expected in the blood vessel or using a fixed number of points (e,g., 
1000). 

10 BRIEF DESCRIPTION OF THE FIGURES 

Non-limiting embodiments of the invention will be described with reference to the 
following description of exemplary embodiments, in conjunction with the figures. The figures 
are generally not shown to scale and any sizes are only meant to be exemplary and not 
necessarily limiting, In the figures, identical structures, elements or parts that appear in more 

15 than one figure are preferably labeled with a same or similar number in all the figures in which 
they appear, in which: 

Fig. 1 is a schematic illustration of a display showing two blood vessel views; 

Fig. 2 is a schematic illustration of a blood vessel showing a stage of a path finding 
method in accordance with an exemplary embodiment of the invention; 
20 Fig. 3 is a schematic illustration of a blood vessel showing a further stage of a path 

finding method in accordance with an exemplary embodiment of the invention; 

Fig. 4 is a schematic illustration of a blood vessel showing a stage of segmentation, In 
accordance with an exemplary embodiment of the invention; 

Fig. S illustrates an optional step of segmentation cleaning, in accordance with an 
26 exemplary embodiment of the invention; 

Fig. 6 illustrates a centerline determination, in accordance with an exemplary 

■ 

embodiment of the invention; 

Fig. 7 is a flowchart of a method of centerline determination, in accordance with an 
exemplary embodiment of the invention; 
30 Fig. 8 is a flowchart of a method of finding an initial path in a blood vessel, in 

accordance with an exemplary embodiment of the invention; and 

Fig. 9 is a flowchart of a method of segmentation based on an initial path, in 
accordance with an exemplary embodiment of the invention. 

4 
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DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS 
Overview of method description and method 

In the following description, a method of finding a centerline is a blood vessel is 
described, in which, Figs. 2-6 illustrate various states during the application of the method, 
5 Fig. 7 is a general flowchart of die method and Figs. 8-9 are more detailed flowcharts of 
particular portions of the method. 

Fig. 7 shows a high-level flowchart 700 for an exemplary embodiment of the invention. 
At 702, a user enters data points indicating a path of the vessel. At 704, a path inside the blood 
vessel is found, connecting the points. At 706, the path, if not found as a single line, is 
10 optionally corrected. Tn an alternative embodiment, this Is done at the end of the process (after 
714), with some or more of acts 708-714 being performed on partial paths. At 708 the blood 
vessel volume is detected (segmentation). At 710, the segmentation is optionally improve! At 
712, a map of the distances inside the blood vessel (from the vessel walls) is determined. At 
714, a centerline is found from the distance map. 
16 Data Entry 

Referring back to Fig. 1, a user enters a plurality of points to indicate the blood vessel 
of interest to the system. It should be noted that in some embodiments, the system will 
generate a centerline for all blood vessels or will automatically detect the blood vessel of 
interest (e.g., based on a constriction therein or based on a desired diagnosis) so that no user 
20 entry is required in such embodiments. 

In a particular implementation, a user' can mark points on one or more of a 3D view, 
and/or one or more slices or projections. Optionally, a windowing function is used to select 
only CT values in a range of possible blood vessel values, thus simplifying the image data set 
and allowing ptojection of relevant anatomical structures. Optionally, the data entry is 
25 performed on a digital image viewing station. 

In the example shown, a user marks a starting point 110 and an ending point 114. A 
middle point 1 12 is optionally provided at a fork. While only three points are shown, more 
points, for example, 4, 6, 10 or more may be used. However, it is a property of some 
embodiments of the invention that only a small number of points is required for entry by a 
30 user, for example, fewer than 1 0, fewer than 8 or fewer than 5, or even 2. 

While the method described below focuses on a single vessel (e.g., 102), optionally, the 
user can mark a bifurcation, for example, by providing a point in a bifurcation vessel 1 04 and a 
point in a bifurcation vessel 106. The method can show a forked centerline. Optionally, the 

5 
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method is applied twice, once for each bifurcation. Optionally, a user indicates when the points 
belong to bifurcations, rather than a convoluted vessel 

Once the method is executed, a centerline 108 is generated Various known imaging 
methods may then be used, for example, methods of slices along the centerline, methods of 
5 slices peipendicular to the centerline and methods virtual endoscopy using the centerline as a 
path and/or viewpoint 

Optionally, the centerline is not restricted to pass through any of the user entered 
points. Alternatively, it may be restricted to some of the points, for example the starting and 
ending points and/or points indicated by the user as crucial. In some embodiments of the 
10 invention it is desirable to receive the user input for where the blood vessel is, but not force the 
user to be very precise in marking points along the exact centerline, which might be time 
consuming and/or difficult 

initial BafeJadiag 

Fig. 8 is a flowchart 800 of a method of finding an initial path in vessel 102, in 
15 accordance with an exemplary embodiment of the invention. Fig. 2 illustrates a stage during 
such path determination. For clarify in this and other figures, only user points 1 10 and 1 12 are 
shorn In general, this path finding method propagates a contour 204 from point 1 10 until it 
meets (e.g., at a point 202) a contour 206 propagated from point 1 12 (or point 1 12 itself). This 
propagation is optionally carried out under a restriction which attempts to reduce the chance of 
20 a point 200 along such a contour from leaving the boundaries of blood vessel 102. Another 
optional restriction is that the path try to have a low cost (e.g., stay near center and/or be 
relatively short). 

A targeted marching algorithm forms the basis of the method of Fig. 8, in accordance 
with an exemplary embodiment of the Invention. In this algorithm, a plurality of points are 
26 stored in a heap, each with an associated cost, and marked as "trial". At each step, a point with 

9 

lowest associated cost is removed from the heap, marked as "alive" and propagated to all its 
neighbors, for example its 6 nearest (Cartesian) neighbors. For each neighbor a cost function 
combining local cost and expected cost to a target are calculated and then the neighbor is put 
in the heap. This is repeated until all the points in the heap are used or the target reached. In 
30 the method of Fig. 8, meeting of points is tracked and if points from two source (user entered) 
points meet, then those points are deemed connected. Points that have not been visited yet are 
marked "far". The different types of marked points ("trial" alive" and "fat") may be processed 
differently. 

6 
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Refcning to Fig. 8, at 802, a heap is initialized and the user points (e.g v 1 1 0, 1 12, and 
1 14) are put in the heap, each with its own unique label and an associated cost of 0. The labels 
are optionally used for finding a starting user point for each point 202, as will be described 
below. Optionally! the user points are ordered geometrically or assumed to be so ordered. The 

6 user points are marked as "trial" and all the other points in the data set are marked as "far". 

At 804 (which may be performed out of order), two user histograms are initialized. As 
will be described below, these histograms help decide if a point under consideration is 
probably inside or outside of a blood vessel. If a different method is used for such 
determination or if no such determination is made, fewer or more histograms may be used, 

10 possibly even no histograms. In an exemplary embodiment of the invention, the two 
histograms used are: 

(a) an "outside" histogram showing die distribution of CT values for the whole data 
set; and 

(b) an "inside" histogram showing the distribution of CT values for points assumed to 
15 be inside the blood vessel. 

In an exemplary embodiment of the invention, the user points are used to initialize the 
"inside" histogram. Optionally the neighbors of these points are entered as well, or other points 
with a high probability of being inside the blood vessel. It should be noted that in this 
implementation the inclusion of "in" points in the "outside" histogram is ignored, but it need 
20 not bo. 

Optionally, one or both histograms are binned, for example to ranges of 5 or 30 CT 

< 

numbers. Alternatively, each gray value added to a histogram is added as a swath of 30 values. 

In the description the term histogram has been used to stand for various methods of 
tracking data value distributions. However, it is not necessary to store a histogram. For 

25 example, a data vector may be used. Alternatively, a curve which approximates the distribution 
may be used. In another example, a neural network is used. Further, other methods of 
determining likely hood may be used, for example, as is known in the ait 

At 806, a point with a lowest associated cost is drawn from the heap. If it id not a user 
point (or even if it is) it is used to update the "inside" histogram", by adding the point's gray 

30 value as a valid blood vessel value. Optionally, this allows points (hat arc determined to be 
outside the blood vessel (e.g., high cost) to be left out of the histogram. Alternatively, the 
"inside" histogram may be updated as an unlabeled point is reached In one such 
implementation, a value is entered into the "inside" histogram with a weight indicative of the 

7 
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probability that that value is actually inside. For example, a Yaiue Is entered with a weight 
based on the cost function (the probability of being outside the blood vessel) described below. 
Alternatively or additionally, for the first 1.00 (or other number) points (e.g., voxels) 
encountered, the neighboring 50 gray levels are entered for each point, with a weight that 

5 decreases as a function of the number of points encountered so far (e.g.. a linear reduction 
from 100 to 1). These and/or other methods may be used to bias the histogram so that gray 
values that are outside the histogram are prevented from taking oyer the histogram. If 
additional data (e.g., base don the medical procedure or expected diagnosis) regarding what 
values arc expected inside or outside blood vessels is available, h may be used to bias the 

1 0 histogram. Alternatively, it may be applied during the cost calculation described below. 

Generally, a drawn point is marked as "alive" and will not be returned to the heap. 
However, in other implementations, such a point may be returned, for example based on 
finding a lower cost path (and/or estimate) thereto. 

At 808, the neighbors of the point are identified. It should be noted that the data set 

1 5 being used can be a 3D data set of cubic voxels. Each point thus has 26 potential neighbors. In 
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As is known in the art of computing, the neighbors may all be determined at once, or one by 

one, for example as a previous neighbor is process, or even in parallel. For other 

representations (e.g., eight sided voxels), a different number of neighbors maybe considered. 
20 In addition, also for cubic voxels, a different number of neighbors, for example, fewer or 

greater than six, may be considered. 

For each such considered neighbor, references 810-828 are applied, and then a next 

point is drawn from the heap (806 again). A stopping condition may be applied, for example at 

a reference 830, as will be described in more detail below. 
25 Referring to Fig. 2, the point drawn from the heap can be, for example point 202 and its 

neighbors be indicated by a reference line 208. 

At 810, a check is made to see if the neighboring point is "alive", "fei" or "trail". If it is 

"alive", a check (812) whether it has the same label as point 202. If it has the same label (814), 

the neighbor point is ignored. 
30 If the neighbor point did not have the same label as point 202, a note is made that a 

connection between two user points was made (816). Once a uBer point is connected on either 

side to other user points (or, if an end point, on only one side), propagation from that point is 



8 
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optionally stopped Once all user points are connected, the process can be stopped (e.g., at 
830). to either case, a substantial reduction in the task is optionally achieved (818). 

If the neighbor point is "far", it is labeled (at 820) with the label of point 202, and 
marled as "trial". In some implementations! a same flag is used both for the label and for the 

6 marking of a pint and "far", "alive" or "trial". In some implementations, it is the presence in 
the heap that marks a point as "trial" and the distinction between "alive" and "far" is whether 
or not the point has a label. 

At 822, a local cost of the neighbor point is calculated, for example, using the formula 
cost - p(value|out)/(p(value|in)f p(value|out)), where p(value|in) answers the question "what is 

10 the chance of seeing this value given that we are inside of the vessel" The p0 function is 
optionally based on the "inside" and "outside "histograms". This is a private case of Baystan 
probability calculation in which the weighing is even fOT being inside or outside. In some 
embodiments, non-even weighting is sued. Optionally, prior knowledge about the data set may 
be used to skew the weighting or to change the pQ functions. This local cost indicates the 

15 likelihood of the point being outside the blood vessel. 

It should be noted that point 200, for example, being at the boundary of the blood 
vessel will probably have neighbors with very high costs (e.g., as they are outside the blood 
vessel) and thus they will probably not be drawn. Also, it may be expected in some situations 
that as a point is closer to the wall of the blood vessel it will be more likely to have a value 

20 associated with being outside die blood vessel. However, this is not always correct (e.g., due to 
non-uniform mixing of a contrast material). 

In an exemplary embodiment of the invention, the resulting local cost for the neighbor 
point is a weighted average of the cost formula for points in its vicinity. In an exemplary 
embodiment of the invention, the cost formula is applied to all 26 nearest neighbors and to the 

25 point itself and averaged using a weighing as follows: 4 for center point, 3 for Cartesian 
neighbors, 2 for semi-diagonal and 1 for diagonals. Then the result is normalized to the range 
0..1 (e.g. f by dividing by 54). In alternative embodiments of the invention, fewer or more 
immediate (or less immediate) neighbors may be averaged 

At 823, the local cost is used to calculate a path cost. In an exemplary embodiment of 

30 the invention, the method used in fast marching is used also for targeted marching. In one 
example, the following equation is solved for u, a searched for path cost 

9 



« 

■ 

032/03892 

(max(«-<7 w _ 1 .«-l/ w+1 ,o) 

where P is the local cost and U are the path costs for each of the points that are 
Cartesian neighbors. In some cases a solution is not possible. A best fit may be searched for. 
Alternatively, one or mote of the "max" units may be replaced by zero. In an exemplary 

5 embodiment of the invention, a "max" unit to be replaced by zero is selected in the following 
manner. For each "max" unit, the smaller U value is found. Then, the "max" unit for which the 
smaller U value is largest, is selected for removal. 

Points for which no path cost has been calculated, are optionally assumed to have 
infinite cost. While the description suggests storing the total cost, a path cost may also be 

10 stored for each point Alternatively, the path cost may be extracted from the total cost 

At 824, a target point from the user points is selected, for which to calculate the 
estimated future cost In some cases a .target is available on only one side and a nearest one is 
selected. In an exemplary embodiment of the invention, the following method is used to select 
between two points. A vector v connecting the neighbor point and the user point it started from 

15 is projected onto two normalized vectors tl and t2 which connect the starting point with the 
candidate target points. The point whose vector has the larger projection, is selected as a target 
point This can be calculated by finding which of <v,tl> and <v,t2> is larger (where <,> 
denotes a scalar multiplication of vectors). 

At 826, a cost to the target is estimated. In an exemplary embodiment of the invention, 

20 the estimate is based on the distance to the target and the average cost per unit distance so far. 
In one implementation, the average unit cost is the path cost to the neighbor point divided by 
the distance from the starting user point to the neighbor point. Then, the estimation of the 
future cost is optionally generated by multiplying the average unit cost by the geometrical 
distance to the target point 

26 At 828, the path cost and the estimated cost are added together to give a total cost, 

which is associated with the point and added with the point to the heap. For neighbor points 
which arc trial points (e.g., in the heap), the associated path cost is updated (if lower), as the 
point docs not need to be added. 

10 
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Optionally, the cost of a point is normalized to the distance between the starting user 

point and the target user point 

It should be noted that several types of costs are provided in some embodiments of the 

invention: 

a) Path cost - the cost to reach a point 

b) Local cost - the cost added by particular properties of the point, for example CT 
number, 

c) Estimated cost - the cost estimated for reaching a target point. 





m 


1 


hi 



10 associated with each point 

At 830, various stopping conditions may be tested. For example, if the lowest cost 
point in the heap is too high or if too much time has passed the path finding method may be 
stopped. In an exemplary embodiment of the invention, a budget of points to be visited is 
provided for various stages of the method, for example, 2,000,000 points each for path finding 

15 and segmentation. 

At 832, the initial path is generated Fig. 3 shows a vessel segment in which a plurality 

of contours 204 were propagated and showing an initial path 300 connecting user points 110 

and 112. This line may be found, for example, by connecting together all the meeting locations 

between different labeled points using, for example, (generally less preferably) straight lines or 
20 lines that travel along the lowest cost path (e.g., using gradient descent or other methods). In 

an exemplary embodiment of the invention, this path is found by advancing the path always to 

the neighbor with the lowest cost 

It should be noted that many path finding methods may be used instead of the method 

of flowchart 800. In an exemplary embodiment of the invention, however, a method is chosen 
25 in which any such path remains inside the blood vessel volume. Desirably, for reducing 

processing costs and/or preventing the determination of bad paths, the search is also focused to 

remain inside the blood vessel 

Path correction 

At 706 (Fig. 7), the found initial path may (optionally) be corrected. In one instance, a 
30 user may desire to correct the path. In another instance, a fall path may have not been found 
(e.g., due to a time limit or the existence of an occlusion in the blood vessel). In an exemplary 
embodiment of the invention, a user can connect between found partial paths. 
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In an exemplary embodiment of the invention, partial paths are connected by applying 
(he targeted marching method with a local cost of I. This will tend to connect segments with 
straight lines. The points to be attached, can be, for example user points, or points on the 
contours extended from nser points 

Alternatively or additionally, the targeted marching method may be reapplied for 
sections where a path was not found, using relaxed constraints, for example, accepting more 
gray values as possibly being part of blood vessels. 

Alternatively, a plurality of partial paths resulting ftom the method of Fig. 8, are passed 
on to the following stages and the final generated centerline is composed by linking together 
centedines generated for each segmentation corresponding to a partial path. 
Segmentation 



30 









111 





data set is the blood vessel of interest In effect, this selection defines a segment that is the 
blood vessel (or a section thereof between the user points 1 10 and 1 14). 

Fig, 4 shows a section of blood vessel 102, shown radially curved, for clarity, during 
such segmentation. Fig, 9 is a flowchart 900 of a method of segmentation based on initial path 
300, in accordance with an exemplary embodiment of the invention. In an exemplary 
embodiment of the invention, this method expands initial path 300 to fill the blood vessel 
volume. Flowchart 900 is of a method based on a "fast marching" method, for example as 
described in Laurent D. Cohen, R. Kimmel, "Global Minimum for Active Contour Models: A 
Minimal Path approach", in International Journal of Computer Vision 24(1) (1997), 57-78 and 
in J. A, Sethian, "Level Set and Fast Marching Methods**, Cambridge university press (1996), 
the disclosures of which is incorporated herein by reference. However, other methods can be 
used as well, for example morphological methods and expansion of centerline methods. In 
addition, in an exemplary embodiment of the invention, various methods to prevent or reduce 
"leaking" from the blood vessel are optionally practiced, for example as described below. It 
should be noted that the fast marching method does not include target selection and target cost 
estimation, but may use a same type of conversion form local cost to path cost 

At 902, the heap for the fast marching method is initialized, with all of the points along 
path 300, optionally, with their immediate neighbors, all with a cost of 0. 

At 904, each point is optionally associated with a parameter. This may be used to 
prevent leaking and/or to ensure correct direction of propagation, as described below. In an 
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exemplary embodiment of the invention, the parameter is an integer ordinal number starting at 
1 for point 110. 

At 906, the blood vessel/initial path is partitioned into a plurality of portions (shown as 
406, 408, 410 and 412). In an exemplary embodiment of the invention, each portion is made 
5 SO points wide along path 300. Optionally, the size of the segments is a parameter which may 
be related, for example to the expected vessel diameter. 

At 908, a local "inside" histogram is built for each portion 406-412. Optionally, the 
local histograms are smoothed and/or shared between nearby portions, for example, a same 
local histogram may be shared between three segments. Optionally, the histograms are binned, 
10 Alternatively, each gray value of path 300 is inserted together with values 20 higher and 20 
lower (or some other range). Alternatively, no local histogram is used, for example, reusing the 
"inside" histogram of Fig. 8 or not using a histogram at all. 

At 910, a plurality of boundary planes 400, 402 and 404, which separate the portions, 
are defined. In an exemplajy embodiment of the invention, these planes are generally 
15 perpendicular to the path 300, and as a result somewhat perpendicular to the vessel walls, 
however, that is not guaranteed In an exemplary embodiment of the invention, the following 
definition is used. Vi is a vector that connects the two ends of a path section (normalized). Ui, 
is a vector that defines the boundary plane between a portion i and a portion i+I, and is an 
average of Vi and Vi+1. The boundary planes are petpendicular to Ul. Ui are optionally 
20 normalized. 

■ 

At 912, a point with a lowest associated cost is drawn from the heap. Optionally, each 
portion has a separate heap, but this is not required 

At 914 rubrics 916-923 are applied to the neighbors of the drawn point Optionally, 

■ 

only the six Cartesian (nearest) neighbors are considered. In an exemplary embodiment of the 
25 invention, "alive" neighbor points are ignored, while "trial" neighbor points are allowed to 
have their cost updated in the heap. 

For clarity, a neighbor point 418 is discussed below. 
At 916, a local cost for point 41 8 is calculated (described in greater detail below). 
At 918, the value of point 418 is optionally used to update the local histogram. The 
30 methods described above are optionally used,- with a smaller cut-off (e.g., special treatment for 
the first 20 points). Alternatively, the histograms need not be updated, just being based on the 
original path 300 and its neighbors (e.g„ nearest 6 or 26 fox different embodiments). 
At 920, a parameterization value is associated with point 418, 

13 
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At 922 a decision is made if to put point 41 8 into the heap or not 

At 923, a point may be added to the heap or its cost updated, if so decided at 922. 

Rubrics 912-923 are repeated. Optionally, a stopping decision 924 may be applied. One 
example stopping decision is to stop if too many steps were performed (e.g., too many points 
5 considered). Another example is a time limit Another example is if the heap is empty, 
Optionally, the stopping conditions are applied separately for different path portions. 

Referring back to 916, the cost function optionally includes a logic portion and a 
calculated portion. The calculated portion decides a value for the cost and the logic portion 
decides if to give a binary cost value (e.g., 0 or infinite). In an exemplary embodiment of the 
10 invention, the cost function is as follows: 

(a) Determine the path portion point 418 belongs to, for example, based on its 
parameterization value, or the parameterization value of the point it originated from. 

(b) Collect and count the number of nearby points that are semi-diagonal to it These 
points are neighbors that share exactly one of an X, Y and Z plane used to locate point 418 in 

16 space. P is the average local cost (e.g., calculated using the formula 
cost=p(valuelout)/(p(value|in)^(valuelout)) t based on the local histogram) for each of the 
collected points. 

(c) If AD, the number of collected points that are "alive" (e.g., already processed) is 
smaller than 2, the local cost for point 418 is infinite. Optionally, this prevents leakage and/or 

20 smoothes the resulting segmentation. 

(d) If AD is greater than 6, then the local cost for point 418 is 0. 

(e) If the average local histogram cost P is greater than 0.6, then the local cost is 
infinite. Optionally, this helps avoid leakage from the blood vessel. 

(f) Otherwise, the local cost is a function of P and AD, that increases with P and 
26 decreases with AD. In an exemplary embodiment of the invention, P is replaced by a weighted 

average of the local cost for point 418 (20%) and the local costs for die semi-diagonals (80%). 

It should be noted that this particular cost function generally prevents propagation with 
high curvature, as increased AD is associated with small local curvature, This means that 
segmentation method prefers filling cavities over generating spikes. 
30 At 917, a path cost is calculated from a local cost, for example using a fast marching 

method, for example, as described above in 823. 

Referring back to 920, a parameter value for point 418 is calculated. In an exemplary 
embodiment of the invention, the new value is selected so that the gradient of parameterization 

14 
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of the points in general is generally parallel to the gradient of parameterization of path 300, In 
. an exemplary embodiment of the invention, this is achieved by selecting a value such that the 
scalar multiplication of a parameterization gradient and a cost gradient is zero. Both gradients 
are optionally interpolated or extrapolated from the points used to calculate the local cost (e.g., 
5 the AD points). Alternatively, a different set of vicinity points is used, for example, the six 
nearest neighbors or all the immediate neighbors. 

Referring back to 922. A decision not to put a point in the heap is optionally based on 
one or more of the following conditions: 

(a) If the cost is infinite or above some predefined threshold. 
10 (b) If point 41 8 goes beyond the planes at the path portion ends. Optionally, a safe zone 

(e.g., 414) is defined within which points from one path portion can encroach on another path 
portion. In an exemplary embodiment of the invention, the origin portion of point 418 is 
determined from its parameterization value. In an exemplary embodiment of the invention, the 
encroachment is determined as follows. A vector connecting point 418 to each of the ends of 
1 5 the portion is projected on each of the plane normal vectors U of these ends. If the projections 
are positive (or above a threshold, for example -20, to accommodate zone 414), then point 418 
is within its allowed portion. 

(c) If point 418 strays too far from path 300. This may serve to prevent leakage or 
travel into side branches of blood vessel 300, for example a side branch 420. In general, it is 
20 assumed that the width of the blood vessel changes gradually. Thus, any sharp changes in the 
profile indicate a leakage or excursion into a side vessel. While this may be dealt with during 
post processing, In an exemplary embodiment of the invention, it is dealt with as part of the 
segmentation. 

In an exemplary embodiment of the invention, a count is maintained of the number of 
25 points with each parameterization value (that is points with the same integer portion of a 
parameterization value). Any points with a count that is 150% of an average count avcount, are 
discarded. Any points with a count between' 100% and 150% are penalized by multiplying 
their cost by count/avcount Optionally, the average count avcount is made on the 9 parameter 
values greater and 9 parameter values smaller than that of point 41 8. 
30 Alternatively or additionally, for example when segmenting blood vessels in the neck, 

a similar method is applied to the path costs. Points with a path cost that is greater than 1 50°/o 
the average path cost of points with the same parameterization value (or an average of 9 points 
on either side) are rejected from the heap. And points between 100% and 150% are penalized. 

15 
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(d) If point 418 was not propagated generally perpendicular to path 300. For example, 
assuming point 418 has a same parameterization value (integer portion) as a point 416 on path 
300, the angle formed between point 418, point 416 and point 110 should be close to 90 
degrees. Optionally, this is applied in neck cases. 

5 In an exemplary embodiment of the invention, the propagation condition is estimated 

as being the scalar product of a vector connecting points 416 and 418 and a gradient of the 
path cost If the product value is above 0.5, the angle is considered to be too great 

The following features of flowchart 900 should be noted as being possibly useful in 
preventing leaks: avoidance of areas with improbably gray values; avoiding narrow paths by 

10 controlling local curvature; avoid undue propagation along the path using boundary plane; 
limiting sideways leakage by forcing smoothness and/or limiting the total number of 
propagation steps. Optionally (e.g., for neck images) an infinite cost value is provided for 
points with an unlikely gray value. For example, if the likelihood of being outside is over a 
threshold such as 90%, the cost is made infinite. 

15 Clean segmentation 

Optionally, at 710, the segmentation is cleaned. In an exemplary embodiment of die 
invention, small leaks are corrected by removing all points that have two neighbors on 
opposing directions that are outside the segmentation. This is shown in Fig. S, where a small 
leak 500 is removed using this method. 

20 Distance transform 

Once the segmentation is found, a centerline may be determined using one or more of 
various methods. In an exemplary embodiment of the invention, the centerline is determined 
by first generating a map of the distances from the boundaries of the segmentation and then 
finding a path along the distances. Alternatively, other methods, such as morphological 

25 skeletonization or thinning arc used. 

In an exemplary embodiment of the invention, a fast marching method is used to 
associate with each point in the segmentation a distance from the outside of the segmentation. 
The fast matching method is applied with each point in the boundary of the segmentation 
being initially in the heap with a cost of 0 and the local cost function being 1 (e.g., a pure 

30 geometrical path cost calculation). Other methods may be used as alternatives, possibly 
inferior, to fast marching, for example, Dijkstra's method of breadth-first search. In an 
exemplary embodiment of the invention, points on the boundary of the segmentation are 

16 



Copy provided bv USPTO from the 1FW Imaqe Database on 01/12/2005 



032/03892 

defined as any points with neighbors outside the segmentation. Fig. 6 shows a plurality of lines 
600 each corresponding to an iso-cost set of points in the data set 
Centerline finding 

Referring to Fig. 6, various methods may be used to find a path from point 1 10 to point 
6 1 12. For example, a ridge following method may travel along the "ridge" of highest distance 
values found in 712, optionally with some smoothing. 

In an exemplary embodiment of the invention, a targeted marching method is used, 
such as described above for initial path finding. In this implementation, only the two end 
points (110 and 114) are placed in the heap, with the target at any step being the end point 
10 from which the method did not start. While even a single point could be used, a potential 
advantage of using two points is that search time may be reduced. 

While the local cost function can simply be the inverse of the distance value found in 
712, in an exemplary embodiment of the invention, a more complex function is used, which 
may assist in providing a smooth centerline. Optionally, this Amotion takes into account the 
16 diameter of the vessel optionally enabling the method to tradeoff a desire to stay in a center of 
the blood vessel with a desire to reduce curvature. 

In an exemplary embodiment of the invention, the Mowing method is used to 
calculate the local cost function of a point P: 

(a) p is set to be the parameterization value found in 708. 
20 (b) MaxD(P) is an estimate of the diameter at point p, for example, the maximum 

distance achieved for any point with parameter p. 

(c) avMaxD(P) is an average on the neighbors of p, for example, 3 higher and three 

lower. 

(d) D(v) is die distance of point P, found in 711 
26 (e) cost= alpha*beta A (D(P)/avMaxD(P))+omega 

(f) omega is a small value used for smoothing, generally close to zero. Alpha is a 
scaling term selected to adjust the other values (e.g., to match omega) and beta is a value 
reflecting a tradeoff between desiring a short path and desiring the path go in a middle of the 
blood vessel. Generally, shorter paths have a greater curvature. Optionally, normalizing D(P) 
30 makes the tradeoff of center on curvature more independent of the blood vessel diameter. In 

* 

one example, omega is 0.001, alpha is 10,000 and beta is 0.00001. 

In some embodiments of the invention, a single path was not found in at 704, but rather 
a plurality of partial paths. Each such partial path is processed as described above and at the 

17 
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end, their centcrlines 600 connected, for example, using straight lines. Optionally, smoothing 
is provided at the connection points. Optionally, such smoothing is applied to the attachment 
section, not to the calculated centerline. 

The above has been described for CT images. However, it may be applied to other 
S modalities as weD, albeit that such modalities may not be as problematic. In such modalities 
various of the parameters may need to be changed and/or optimized. For example, the degree 
of smearing in histogram values may be changed. 

The above method may be applied on vessels of various sizes, for example from 2 mm 
to 25 mm, for example with a voxel size of less than 1 ram or less than 0.7mm or 0.5 mm. For 
10 example, a centerline may be found for the aorta, from its arch to the iliac arteries, for renal 
arteries, for vertebral arteries, for internal carotids. 

While the method may be applied to voxel based methods, other representation 
methods can also be used with the vessel. 

It should be appreciated that while a complete method of finding a centerline i$ 

15 described, in some exemplary embodiments of the invention, parts of the above method are 
used for other image processing task, possibly on other body parts and possibly on non- 
medical images. However, various of the above described features appear to be capable of 
providing a beneficial result for the task of centerline finding in blood vessels. 

It should be noted that not all of the problems mentioned in the background can be 

20 guaranteed to be solved by all embodiments of the invention, instead, suspect situations may 
be provided to a physician to decide on, for example, indicating areas where a centerline could 
not be found Alternatively or additionally, a physioian can add points and/or change various 
parameters and ask the method to be executed again. 

It should be noted that in some embodiments of the invention, not all of the user 

25 entered points are on the centerline. For ©cample, possibly only the two end user points arc 
used. Even these points are not required to be on the centerline, in some implementations. In 
an extreme case, a U9er can indicate one point in a blood vessel and allow the method to extend 
its search to the edge(s) of the volume. In another extreme case, all the blood vessels in a 
volume are centerlined, for example, based on an otherwise identified blood vessel. 

30 In the above description many numbers are provided. In some implementations these 

numbers are parameters. In general, a user can optimize the values and/or other properties of 
the algorithm for particular cases, for example, particular types of CT images. In any case, this 
numbers should be considered examples of possible numbers and not absolutely limiting on 
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the invention. In particular, some parameters may be changed for particular ranges of blood 
vessel diameters and/or types of problems. Also, when transferring to other problems, such as 
virtual endoscopy (or viewing of other lumens or solid elongate organs), other values may be 
useful. 

5 A particular property of some embodiments of the invention is that only a relatively 

small number of points is required, and that other than requiring these points to be inside the 
vessel, they may be loosely placed (in some implementations, the two end points are held to 
higher standards). In an exemplary embodiment of the invention, a blood vessel requires fewer 
than 2+3N points, where N is the number of occlusions and/or number of major stenosed areas 
10 in a blood vessel. In some cases, only 2+2N points are required. In others, 2+N points are 
sufficient 

It should be noted that various types of vicinities and neighborhoods were used above. 
In other implementations, the number of points in certain neighborhood may be increased or 
reduced. In one example, where heavy processing is required, fewer neighbors are considered. 
16 In another, example, for example for traveling, either 6 or 26 neighbors are used, to prevent 
problems caused by semi-diagonals and their neighbors. It should alsu be noted that while 
averaging is widely practiced, it not required in every case, while it does usually provide a 
generally desired damping effect 

The present invention has been described using non-limiting detailed descriptions of 
20 embodiments thereof that are provided by way of example and are not intended to limit the 
scope of the invention. It should be understood that features and/or steps described with 
respect to one embodiment may be used with other embodiments and that not all embodiments 
of the invention have all of the features and/or steps shown in a particular figure or described 
with respect to one of the embodiments. It is noted that some of the above described 
25 embodiments may describe the best mode contemplated by the inventors and therefore include 
structure, acts or details of structures and acts that may not be essential to the invention and 
which are described as examples. 

Structure and acts described herein are replaceable by equivalents which perform the 
same function, even if the structure or acts are different, as known in the art. Hierefore, the 
30 scope of the invention is limited only by the elements and limitations as used in the claims. 
When used in the following claims, the terms "comprise", "include", "have" and their 
conjugates mean "including but not limited to". 

> 
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While the application has been described as methods! it is meant to also encompass 
apparatus for carrying out the invention, for example, suitably programmed servers and client 
computers and/or media having suitable software thereon, for example, diskettes and computer 
RAM. 
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CLAIMS 

■ 

1. A method of centerllne determination for a blood vessel in a data set, comprising: 
5 receiving at least one start point and one end point inside a blood vessel volume; 

automatically determining a path between said points that remains inside said volume; 
automatically segmenting said blood vessel using said path; and 
automatically determining a centerline from said segmentation. 

10 2. A method according to claim 1, wherein automatically determining a path comprises 
determining using targeted marching. 

• - • 

3. A method according to claim 2, wherein targeted marching uses a cost Amotion 
incorporating both path cost and estimated future path cost 

15 

4. A method according to claim l v wherein said automatically segmenting comprises 
segmenting while reduces leakage out of said volume during said segmenting. 
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